Using Geoscience Australia WOfS algorithm.
This notebook requires the "WOfS Environment", which needs to be installed on the first run of this notebook by executing the following cell. Once done, set the environment in the top-right corner.
# Install the WOfS environment - Only needed on the very first execution of this notebook
!../tools/install_wofs.sh
WOfS Environment kernel already installed
product = "landsat8_c2l2_sr"
longitude = (119.8242517, 120.0350519)
latitude = (-4.2013799, -3.9445384)
time = ('2020-01-01', '2020-12-31')
output_crs = "EPSG:32650"
resolution = (30, -30)
# Where to save the DEM fetched in ODC
DEM_PATH = "/home/jovyan/dems/srtm_lake_tempe.tif"
from os import environ
environ["AWS_HTTPS"] = "NO"
environ["GDAL_HTTP_PROXY"] = "easi-caching-proxy.caching-proxy:80"
print(f'Will use caching proxy at: {environ.get("GDAL_HTTP_PROXY")}')
Will use caching proxy at: easi-caching-proxy.caching-proxy:80
import sys
sys.path.append("../../hub-notebooks/scripts")
from pathlib import Path
import xarray as xr
import rioxarray
from app_utils import display_map
from wofs.virtualproduct import WOfSClassifier
from datacube import Datacube
from datacube.utils.aws import configure_s3_access
configure_s3_access(
aws_unsigned=False,
requester_pays=True,
);
dc = Datacube()
# Display available products
products_info = dc.list_products()
products_info
| name | description | license | default_crs | default_resolution | |
|---|---|---|---|---|---|
| name | |||||
| copernicus_dem_30 | copernicus_dem_30 | Copernicus 30m Digital Elevation Model (GLO-30) | None | None | None |
| landsat5_c2l2_sr | landsat5_c2l2_sr | Landsat 5 Collection 2 Level-2 Surface Reflect... | None | None | None |
| landsat5_c2l2_st | landsat5_c2l2_st | Landsat 5 Collection 2 Level-2 UTM Surface Tem... | CC-BY-4.0 | None | None |
| landsat7_c2l2_sr | landsat7_c2l2_sr | Landsat 7 USGS Collection 2 Surface Reflectanc... | None | None | None |
| landsat7_c2l2_st | landsat7_c2l2_st | Landsat 7 Collection 2 Level-2 UTM Surface Tem... | CC-BY-4.0 | None | None |
| landsat8_c2l1 | landsat8_c2l1 | Landsat 8 Collection 2 Level-1 (top of atmosph... | CC-BY-4.0 | None | None |
| landsat8_c2l2_sr | landsat8_c2l2_sr | Landsat 8 Collection 2 Surface Reflectance, pr... | None | None | None |
| landsat8_c2l2_st | landsat8_c2l2_st | Landsat 8 Collection 2 Level-2 UTM Surface Tem... | CC-BY-4.0 | None | None |
| landsat9_c2l1 | landsat9_c2l1 | Landsat 9 Collection 2 Level-1 (top of atmosph... | CC-BY-4.0 | None | None |
| landsat9_c2l2_sr | landsat9_c2l2_sr | Landsat 9 Collection 2 Surface Reflectance, pr... | CC-BY-4.0 | None | None |
| landsat9_c2l2_st | landsat9_c2l2_st | Landsat 9 Collection 2 Level-2 UTM Surface Tem... | CC-BY-4.0 | None | None |
| lpdaac_nasadem | lpdaac_nasadem | NASADEM Merged DEM Global 1 arc second V001 | None | None | None |
| nasa_aqua_l2_oc | nasa_aqua_l2_oc | NASA MODIS-Aqua L2 Ocean Color, regridded to W... | None | EPSG:4326 | (0.01, 0.01) |
| nasa_aqua_l2_sst | nasa_aqua_l2_sst | NASA MODIS-Aqua L2 Sea Surface Temperature, re... | None | EPSG:4326 | (0.01, 0.01) |
| s2_l2a | s2_l2a | Sentinel-2a and Sentinel-2b imagery, processed... | None | None | None |
display_map(x=longitude, y=latitude)
We load the SRTM DEM data from ODC using the lpdaac driver. The Lake Tempe area is mostly flat so the DEM can be made optional.
# Ignore warnings in output
import warnings
from sqlalchemy.exc import SAWarning
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=SAWarning)
from os import environ
from cartopy.crs import PlateCarree
from datashader import reductions
import holoviews as hv
import hvplot.xarray
import matplotlib.pyplot as plt
dem = dc.load(
product="lpdaac_nasadem",
latitude=latitude,
longitude=longitude,
output_crs="epsg:4326",
resolution=(-1/3600, 1/3600),
)
elevation = dem.elevation
options = {
'title': 'Elevation',
'width': 800,
'height': 500,
'aspect': 'equal',
'cmap': plt.cm.terrain,
'clim': (0, elevation.max().values.item()), # Limit the color range depending on the layer_name
'colorbar': True,
'tools': ['hover'],
}
plot_crs = PlateCarree()
elevation.hvplot.image(
x = 'longitude', y = 'latitude', # Dataset x,y dimension names
crs = plot_crs,
rasterize = True, # If False, data will not be reduced. This is slow to load but all data is loaded.
aggregator = reductions.mean(), # Datashader calculates the mean value for reductions (also first, min, max, las, std, mode)
precompute = True, # Datashader precomputes what it can
).opts(**options).hist(bin_range = options['clim'])
dem_path = Path(DEM_PATH)
dem_path.parent.mkdir(parents=True, exist_ok=True)
elevation.rio.to_raster(dem_path)
Load the data from ODC and rename bands as needed by the WOfS classifier.
# Ignore SAWarning in output
measurements = ['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'pixel_qa']
data = dc.load(
product=product,
longitude=longitude,
latitude=latitude,
time=time,
output_crs=output_crs,
resolution=resolution,
measurements=measurements,
dask_chunks={'time': 1},
)
data = data.rename({
"blue": "nbart_blue",
"green": "nbart_green",
"red": "nbart_red",
"nir": "nbart_nir",
"swir1": "nbart_swir_1",
"swir2": "nbart_swir_2",
"pixel_qa": "fmask",
})
data
<xarray.Dataset>
Dimensions: (time: 19, y: 951, x: 785)
Coordinates:
* time (time) datetime64[ns] 2020-01-09T02:10:26.574371 ... 2020-1...
* y (y) float64 -4.65e+05 -4.65e+05 ... -4.366e+05 -4.365e+05
* x (x) float64 8.371e+05 8.37e+05 ... 8.136e+05 8.136e+05
spatial_ref int32 32650
Data variables:
nbart_blue (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
nbart_green (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
nbart_red (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
nbart_nir (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
nbart_swir_1 (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
nbart_swir_2 (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
fmask (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
Attributes:
crs: EPSG:32650
grid_mapping: spatial_reftransform = WOfSClassifier(c2_scaling=True, dsm_path=DEM_PATH)
# Compute the WOFS layer
wofl = transform.compute(data)
wofl
<xarray.Dataset>
Dimensions: (y: 951, x: 785, time: 19)
Coordinates:
* y (y) float64 -4.65e+05 -4.65e+05 ... -4.366e+05 -4.365e+05
* x (x) float64 8.371e+05 8.37e+05 8.37e+05 ... 8.136e+05 8.136e+05
* time (time) datetime64[ns] 2020-01-09T02:10:26.574371 ... 2020-12...
spatial_ref int32 32650
Data variables:
water (time, y, x) uint8 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
Attributes:
crs: EPSG:32650# Uncomment the following line to display WOfS for each time slice. This may take a few minutes
# wofl.water.plot(col="time", col_wrap=5);
The Water Observations Summaries based on https://github.com/opendatacube/odc-stats/blob/develop/odc/stats/plugins/wofs.py are made up of:
count_clear: a count of every time a pixel was observed (not obscured by terrain or clouds)count_wet: a count of every time a pixel was observed and wetfrequency: what fraction of time (wet/clear) was the pixel wet# Rename dimensions as required
wofl = wofl.rename({"x": "longitude", "y": "latitude"})
from odc.algo import safe_div, apply_numexpr, keep_good_only
wofl["bad"] = (wofl.water & 0b0111_1110) > 0
wofl["some"] = apply_numexpr("((water<<30)>>30)==0", wofl, name="some")
wofl["dry"] = wofl.water == 0
wofl["wet"] = wofl.water == 128
wofl = wofl.drop_vars("water")
for dv in wofl.data_vars.values():
dv.attrs.pop("nodata", None)
# Ignore warnings triggered by time slices without data at all
warnings.filterwarnings("ignore", message="divide by zero encountered in true_divide")
warnings.filterwarnings("ignore", message="invalid value encountered in true_divide")
wofl.wet.plot(col="time", col_wrap=5);
# Helper frunction from https://github.com/opendatacube/odc-stats/blob/develop/odc/stats/plugins/wofs.py
def reduce(xx: xr.Dataset) -> xr.Dataset:
nodata = -999
count_some = xx.some.sum(axis=0, dtype="int16")
count_wet = xx.wet.sum(axis=0, dtype="int16")
count_dry = xx.dry.sum(axis=0, dtype="int16")
count_clear = count_wet + count_dry
frequency = safe_div(count_wet, count_clear, dtype="float32")
count_wet.attrs["nodata"] = nodata
count_clear.attrs["nodata"] = nodata
is_ok = count_some > 0
count_wet = keep_good_only(count_wet, is_ok)
count_clear = keep_good_only(count_clear, is_ok)
return xr.Dataset(
dict(
count_wet=count_wet,
count_clear=count_clear,
frequency=frequency,
)
)
summary = reduce(wofl)
A count of every time a pixel was observed and wet.
summary.count_wet.plot(size=10)
plt.title("Lake Tempe – Wet observations counts for 2020");
A count of every time a pixel was observed (not obscured by terrain or clouds).
summary.count_clear.plot(size=10)
plt.title("Lake Tempe – Clear observations counts for 2020");
What fraction of the time was the pixel wet.
summary.frequency.plot(size=10)
plt.title("Lake Tempe – Wet frequency for 2020");
# Save frequencies to high-res image file
summary.frequency.plot(size=20)
plt.title("Lake Tempe – Wet frequency for 2020")
plt.savefig('/home/jovyan/tempe_wofs_2020.png', dpi=600);
